library(tidyverse)
library(plotly)
library(sf)
library(tigris)
library(leaflet)
library(censusapi)
library(zoo)years <- 2017:2020
quarters <- 1:4
types <- c("Electric","Gas")
pge_17_to_20_elec_and_gas <- NULL
for (year in years) {
if (year == 2020) {
quarters <- 1:2
}
for (quarter in quarters) {
for (type in types) {
filename <-
paste0(
"PGE_",
year,
"_Q",
quarter,
"_",
type,
"UsageByZip.csv"
)
temp <- read_csv(filename)
pge_17_to_20_elec_and_gas <- bind_rows(pge_17_to_20_elec_and_gas, temp)
saveRDS(pge_17_to_20_elec_and_gas, "pge_17_to_20_elec_and_gas.rds")
}
}
}ca_counties <- counties("CA", cb=T, progress_bar = F)
bay_county_names <-
c(
"Alameda",
"Contra Costa",
"Marin",
"Napa",
"San Fransisco",
"San Mateo",
"Santa Clara",
"Solano",
"Sonoma"
)
bay_counties <-
ca_counties %>%
filter(NAME %in% bay_county_names)
usa_zips <-
zctas(cb = T, progress_bar = F)
bay_zips <-
usa_zips %>%
st_centroid() %>%
.[bay_counties, ] %>%
st_set_geometry(NULL) %>%
left_join(usa_zips %>% select(GEOID10)) %>%
st_as_sf()pge_final <-
pge_17_to_20_elec_and_gas %>%
filter(CUSTOMERCLASS %in%
c(
"Elec- Residential",
"Elec- Commercial",
"Gas- Residential",
"Gas- Commercial"
),
ZIPCODE %in% bay_zips$ZCTA5CE10) %>%
mutate(
TOTALKBTU = case_when(
!is.na(TOTALKWH) ~ TOTALKWH * 3.412,
!is.na(TOTALTHM) ~ TOTALTHM * 99.976
)
) %>%
group_by(YEAR, MONTH, CUSTOMERCLASS) %>%
summarize(
MONTHLYKBTU =
sum(
TOTALKBTU,
na.rm = T
)
) %>%
mutate(Date = as.yearmon(paste(YEAR, MONTH, sep = "-")))pge_chart <-
pge_final %>%
ggplot() +
geom_line(
mapping = aes(
x= Date,
y= MONTHLYKBTU,
color = CUSTOMERCLASS),
stat = "identity"
) +
labs(
y = "Monthly kBTU",
title = "PG&E Bay Area Monthly Electricity and Gas Usage, 2017-2020",
fill = "Electricity Type"
)
pge_chart %>% ggplotly()Overall, gas consumption followed a similar trend as electricity consumption, but there are a couple of exceptions. In April of 2020 commercial gas consumption was actually 3.3% higher than in April of 2019, and in May of 2020 residential gas use was 14% lower than in May of 2019. The largest difference in residential gas use occurred in April, when consumption was 31% higher in 2020 than in 2019. Moreover, the largest difference in commercial gas use occurred in May, when consumption decreased ~17% from 2019 to 2020. This trend suggests the opposite conclusion from that of electricity consumption; the first full month of the stay-at-home order had a significant impact on residential gas consumption, whereas commercial gas consumption experienced its greatest change in May.
We cannot definitively attribute these trends to COVID-19, but it’s evident that there are changes in consumption and it’s reasonable to assume that the pandemic is, to a certain degree, playing a role in causing these differences.
Lastly, it’s important to recognize the differences in consumption data for September 2017 from the rest of the trend. It appears that there’s an error in the PG&E data, which is leading to the appearance of spikes in energy consumption for that month. If the data were correct, we could expect to see a trend similar to that of the other years.
pge_res_elec_covid <-
pge_17_to_20_elec_and_gas %>%
filter(YEAR %in% (2019:2020), MONTH %in% (4:6),
CUSTOMERCLASS == "Elec- Residential") %>%
mutate(
ZIPCODE = ZIPCODE %>% as.character()
) %>%
group_by(ZIPCODE, YEAR) %>%
summarize(
TOTALKWH =
sum(TOTALKWH, na.rm = T)) %>%
pivot_wider(
names_from = YEAR,
values_from = TOTALKWH
) %>%
rename(
KWH2019 = "2019", KWH2020 = "2020"
) %>%
mutate(
percent_change =
((KWH2020 - KWH2019)/KWH2019) * 100
) %>%
right_join(
bay_zips %>% select(GEOID10),
by = c("ZIPCODE" = "GEOID10")
) %>%
filter(KWH2019 != 0, KWH2020 !=0, !is.na(KWH2019), !is.na(KWH2020)) %>%
st_as_sf() %>%
st_transform(4326)pge_comm_elec_covid <-
pge_17_to_20_elec_and_gas %>%
filter(YEAR %in% (2019:2020), MONTH %in% (4:6),
CUSTOMERCLASS == "Elec- Commercial") %>%
mutate(
ZIPCODE = ZIPCODE %>% as.character()
) %>%
group_by(ZIPCODE, YEAR) %>%
summarize(
TOTALKWH =
sum(TOTALKWH, na.rm = T)
) %>%
pivot_wider(
names_from = YEAR,
values_from = TOTALKWH
) %>%
rename(
KWH2019 = "2019", KWH2020 = "2020"
) %>%
mutate(
percent_change =
((KWH2020 - KWH2019)/KWH2019) * 100
) %>%
right_join(
bay_zips %>% select(GEOID10),
by = c("ZIPCODE" = "GEOID10")
) %>%
filter(KWH2019 != 0, KWH2020 !=0, !is.na(KWH2019), !is.na(KWH2020)) %>%
st_as_sf() %>%
st_transform(4326)res_pal <- colorBin(
palette = "PiYG",
domain =
pge_res_elec_covid$percent_change,
bins = c(-75, -20, -10, 0, 10, 20, 30)
)
leaflet() %>%
addTiles() %>%
addPolygons(
data = pge_res_elec_covid,
fillColor = ~res_pal(percent_change),
color = "white",
opacity = 1,
fillOpacity = 1.5,
weight = 1,
label = ~paste0(
round(percent_change),
" percent change in ",
ZIPCODE
),
highlightOptions = highlightOptions(
weight = 2,
opacity = 1
)
) %>%
addLegend(
data = pge_res_elec_covid,
pal = res_pal,
values = ~percent_change,
title = "Percent Change in <br> April-June Residential <br> kWH, 2019-2020"
) res_pal <- colorBin(
palette = "RdBu",
domain =
pge_comm_elec_covid$percent_change,
bins = c(-90, -70, -50, -30, -10, 0, 25, 50, 125, 200, 300)
)
leaflet() %>%
addTiles() %>%
addPolygons(
data = pge_comm_elec_covid,
fillColor = ~res_pal(percent_change),
color = "white",
opacity = 1.5,
fillOpacity = 2,
weight = 1,
label = ~paste0(
round(percent_change),
" percent change in ",
ZIPCODE
),
highlightOptions = highlightOptions(
weight = 2,
opacity = 1
)
) %>%
addLegend(
data = pge_comm_elec_covid,
pal = res_pal,
values = ~percent_change,
title = "Percent Change in <br> April-June Commercial <br> kWH, 2019-2020"
)